Role of interactions in a closed quenched driven diffusive system

Driven diffusive systems, owing to their occurrences in large number of physical and biological processes, are of great significance. Some of the familiar experiences such as the flocking of birds or fish [1], ant trails [2, 3], traffic flow [2, 3], or biological transport [4, 5], are just a few examples of such systems. These systems fall into non-equilibrium category which is far less understood than the equilibrium counterpart. However, the systems can settle down to a non-equilibrium steady state (NESS) which has been studied over the years to gain deep insights into their properties [6]. One of the most powerful tools in investigating multi-particle non-equilibrium system is a class of models called the totally asymmetric simple exclusion process (TASEP) [7–11, 30, 35]. It was originally proposed as a simple model for the motion of multiple ribosomes along mRNA during protein translation [12]. The model involves hopping of particles from one site to immediate next site on a one-dimensional lattice with a unit rate and obeying hardcore exclusion principle [9, 12, 13]. Over the years, due to its simplicity, different versions of TASEP have been extensively employed in studies of various aspects of biological motors, vehicular traffic, etc, providing important insights into these complex processes [14–19, 30].

Various versions of TASEP models have been thoroughly investigated. With open boundaries and unit hopping rates, the system settles into one of the three phases depending upon entry rate (α) and exit rate (β). These phases are referred to as high density (HD), low density (LD) and maximal current (MC) [9, 36]. A variant of this system incorporated with weak correlations has a topologically similar phase diagram [20]. In a closed system with unit hopping rates, depending upon the number of particles, HD, LD and MC is obtained [3, 13, 21]. With the particles interacting with energy E, the system exhibits a higher value for MC in case of weak interactions [18]. The unit hopping rates, considered for simplicity, generally do not hold true for majority of systems [3, 17, 22, 23]. For instance, a vehicle on a road can move with non uniform speed; it may slow down when it encounters a sharp turn or a speed breaker, etc, or its speed may be altered due to different speed limits for different parts of the road. Also, in simple TASEP, the particles are non interacting, i.e., hopping rate at a site is constant and remains unaffected by the occupancy of neighboring sites. However, in vehicular traffic it is observed that a vehicle slows down in presence of a vehicle immediately ahead of it, and speeds up if a vehicle behind it starts honking [23]. Thus particles are influence by the presence of another particle. Similarly, experimental studies on kinesin motor proteins, which move along microtubules, indicate that these molecular motors interact with each other [24, 31, 32]. A model for kinesin having nearest neighbor particle interactions has been studied wherein the hopping rates are modified due to interactions taking into consideration the fundamental thermodynamic consistency [18, 25, 33]. Many variants of TASEP have been studied wherein the hopping rates are inhomogeneous [16, 19, 26, 27]. A recent study considers a system consisting of non-interacting particles on a closed lattice with quenched hopping rates i.e., hopping rate depends upon site, and shows the dependence of NESS on total number of particles and minimum of the quenched rates. The minimum acts as bottleneck and a shock in form of localized domain wall is observed in MC phase [17].

Motivated by the relevance of both quenched hopping rate and particle–particle interactions, we incorporate them in the simple TASEP to analyze a generalized version. Taking the recent studies into account, our objective is to answer the following questions. Does the simple mean field (SMF) theory, which worked for the quenched system of non-interacting particles [17], give accurate predictions for the proposed model? If no, then what advanced theory can be applied to obtain the dynamics of such a system? Do the qualitative properties of the system change with the inclusion of the interactions? How does the incorporation of interactions affect the shock phases? We proceed to answer these questions and a few more in the forthcoming sections using various mean-field methods.

We consider a model comprised of N interacting particles on a closed lattice with L sites. Particles move unidirectionally on the lattice and obey the exclusion principle i.e., at most one particle can occupy a single site. A particle hops to the next site only if the target site is empty. Unlike the unit hopping rate in simple TASEP, we consider quenched hopping rates characterized by λ_i for ith site. This incorporation takes into account the inhomogeneity of the paths i.e., twists and turns in microtubules, and speed breakers, speed limits, etc, in roads and thus, makes the model more applicable. Additionally, the particles in the system interact with energy E which is associated with the bonds connecting two nearest neighboring particles, where E > 0 and E < 0 corresponds to attraction and repulsion, respectively [18, 20]. When a particle hops to its next site, it can break or form a bond depending on occupancy of neighboring sites that can be interpreted as opposite chemical reactions. For attractive interactions (E > 0), a particle has a tendency to form a bond with the particle ahead of it, thereby increasing its rate of hopping by a factor q(>1); whereas a particle resists the breaking off from bond with the particle behind it, leading to a decrease in its rate of hopping to next site by a factor r(<1). Similar arguments hold for the repulsive interactions (E < 0). The rates q and r are taken in a thermodynamic consistent manner which defines the forming and breaking of particle–particle bond as q = e^θE and r = e^(θ−1)E , respectively, where θ(0 < θ ⩽ 1) allows the measure of the distribution of energy [20]. Throughout our paper, we assume λ_i = λ(i/L) having 0 < λ(i/L) ⩽ 1 where λ(⋅) is spatially piecewise smooth and has a single point global minimum. Depending upon the occupancy of neighboring sites, the hopping rate from site i − 1 to i is defined as follows (see figure 1):

• When site i − 2 and i + 1 are empty (or occupied), the hopping rate is λ_i .
• When site i − 2 is empty and site i + 1 is occupied, the hopping rate is qλ_i with q ≠ 1.
• When site i − 2 is occupied and site i + 1 is empty, the hopping rate is rλ_i with r ≠ 1.

Zoom In Zoom Out Reset image size

In the absence of interactions i.e., E = 0, we have q = r = 1 and thus we recover the TASEP for closed ring with quenched hopping rate [17]. Further, if λ_i 's = 1, we get the simple closed TASEP [13]. We denote the occupation number of the ith site as τ_i which is assigned 1(0) when it is occupied (unoccupied). Then, the master equation for the temporal evolution of average site occupation number, ⟨τ_i ⟩ (0 ⩽ i ⩽ L − 1), is given as

$$\begin{equation}\frac{\mathrm{d}\langle {\tau }_{i}\rangle }{\mathrm{d}{t}^{\prime }}={J}_{i-1,i}-{J}_{i,i+1}.\end{equation} \tag{ 1 }$$

Here, J_i,i+1 is the particle current from site i to i + 1, given as

$$\begin{align}\hfill {J}_{i,i+1}& ={\lambda }_{i+1}\left(\langle {\bar{\tau }}_{i-1}{\tau }_{i}{\bar{\tau }}_{i+1}{\bar{\tau }}_{i+2}\rangle +r\langle {\tau }_{i-1}{\tau }_{i}{\bar{\tau }}_{i+1}{\bar{\tau }}_{i+2}\rangle +q\langle {\bar{\tau }}_{i-1}{\tau }_{i}{\bar{\tau }}_{i+1}{\tau }_{i+2}\rangle +\langle {\tau }_{i-1}{\tau }_{i}{\bar{\tau }}_{i+1}{\tau }_{i+2}\rangle \right),\hfill \end{align} \tag{ 2 }$$

where ${\bar{\tau }}_{i}=1-{\tau }_{i}$ and i is replaced by i(mod L) whenever i ∉ {0, 1, ..., L − 1}. The four point correlators in equation (2) make it intractable in the present form which, in turn, prompts us to look for an approximation to the correlators. We approach this problem with mean field theory. The basic premise is to break the four point correlators into smaller correlators. We begin with the SMF approximation wherein the idea is to ignore all possible correlations between the particles and probability of products is replaced by products of probabilities, i.e.,

$$\begin{equation}\langle {\tau }_{i}{\tau }_{j}{\tau }_{k}{\tau }_{l}\rangle \approx \langle {\tau }_{i}\rangle \langle {\tau }_{j}\rangle \langle {\tau }_{k}\rangle \langle {\tau }_{l}\rangle .\end{equation} \tag{ 3 }$$

By denoting ⟨τ_i ⟩ by ρ_i and under the approximation given by equation (3), equation (2) becomes

$$\begin{align}\hfill {J}_{i,i+1}& ={\lambda }_{i+1}{\rho }_{i}\left(1-{\rho }_{i+1}\right)\left(\left(1-{\rho }_{i-1}\right)\left(1-{\rho }_{i+2}\right)+r{\rho }_{i-1}\left(1-{\rho }_{i+2}\right)\right.\hfill \\ \hfill & \quad \left.+\enspace q\left(1-{\rho }_{i-1}\right){\rho }_{i+2}+{\rho }_{i-1}{\rho }_{i+2}\right).\hfill \end{align} \tag{ 4 }$$

In thermodynamic limit N → ∞, consider a lattice constant = 1/L. The sites are labeled by x = i and x becomes a pseudo continuous variable between 0 and 1. The temporal variable t' is replaced by t = t'/L. The average density, ρ_i ≡ ρ(x), together with the Taylor series expansion (up to first order derivatives) leads to the following continuity equation:

$$\begin{equation}\frac{\partial \rho }{\partial t}+\frac{\partial J}{\partial x}=0,\end{equation} \tag{ 5 }$$

where

$$\begin{equation}J=\lambda \left(x\right)\rho \left(1-\rho \right)\left(1+\left({\rho }^{2}-\rho \right)\left(2-q-r\right)\right)\end{equation} \tag{ 6 }$$

is the current in the system. In steady state, the average density remains unchanged in time i.e., $\frac{\partial \rho }{\partial t}=0$ due to which equation (5) assures that the current in the steady state is constant. In contrast to the case of simple TASEP with interacting particles, here ρ will not be a constant throughout the lattice [25]. Solving equation (6) leads us to the following expressions for density profile:

$$\begin{equation}{\rho }_{{\pm}}\left(x\right)=\begin{cases}\frac{1}{2}\left[1{\pm}\sqrt{1+\frac{2\left(1-\sqrt{\frac{4J\left(q+r-2\right)}{\lambda \left(x\right)}+1}\right)}{\left(q+r-2\right)}}\right],\quad q,r\ne 1\quad \hfill \\ \frac{1}{2}\left[1{\pm}\sqrt{1-\frac{4J}{\lambda \left(x\right)}}\right],\quad \quad \quad \quad \quad \quad \quad \quad \enspace q=r=1\quad \hfill \end{cases}\end{equation} \tag{ 7 }$$

for all x. Clearly, ρ₋(x) $\left({\rho }_{+}\left(x\right)\right)$ is bounded above (below) by 0.5, and depends upon J. This J can be calculated by using the particle number conservation (PNC):

$$\begin{equation}{\int }_{0}^{1}{\rho }_{a}\left(x\right)\mathrm{d}x=n,\quad a=+,-,\end{equation} \tag{ 8 }$$

where $n=\frac{N}{L}$. For feasible values of ρ(x),

$$\begin{equation}J{\leqslant}\begin{cases}\frac{\left(2+q+r\right)\lambda \left(x\right)}{16},\quad q,r\ne 1\quad \hfill \\ \frac{\lambda \left(x\right)}{4},\quad \quad \quad \quad \quad \enspace q=r=1\quad \hfill \end{cases}\end{equation} \tag{ 9 }$$

for all x. Thus, the maximum possible value of particle current is

$$\begin{equation}{J}_{\mathrm{max}}=\begin{cases}\frac{\left(2+q+r\right){\lambda }_{\mathrm{min}}}{16},\quad q,r\ne 1\quad \hfill \\ \frac{{\lambda }_{\mathrm{min}}}{4},\quad \quad \quad \quad \quad \enspace q=r=1,\quad \hfill \end{cases}\end{equation} \tag{ 10 }$$

where λ_min is the global minimum of the λ(x). At J = J_max, ρ₋(x₀) = ρ₊(x₀) where x₀ is the point of global minimum. For the limiting case, as E → 0, the expressions for density profile and MC agrees well with that obtained in reference [17]. Further if λ(x) = 1, the results exactly match with that of the simple TASEP with non-interacting particles.

From figures 2(a) and (b), it can be observed that the expressions for density obtained from SMF produces exactly the same density profiles qualitatively as well as quantitatively for equal strength of both attractive and repulsive interactions. This outcome is expected due to the fact that q + r is an even function of E [see figure 2(c)]. In order to investigate the validity of the above obtained results, we simulate model by employing Monte Carlo simulations (MCS) [37] and apply random sequential update rules. Here, a site is chosen randomly and is updated in accordance with the system dynamics. The lattice length is taken to be L = 500 and the simulations run for 10⁹ timesteps. To ensure the occurrence of the steady-state, we ignore the first 5% of the timesteps and the average density of particles is computed over an interval of 10L.

Zoom In Zoom Out Reset image size

It can be readily observed from figures 2(a) and (b) that, in contrast to SMF, MCS shows that the effect of attractive interaction is distinct from that of repulsive interactions for equal strength. For a fixed value E, the density profiles obtained from SMF differs from MCS for any value of n. For instance, in figures 3(a) and (b) it can be observed that SMF assures the presence of a shock which is not observed in MCS result; whereas in figure 3(c), SMF does not predict a shock in density profile as opposed to MCS. Furthermore, when both SMF and MCS predict same phase for some fixed value of n and E, the quantitative difference in the density profile can be seen as shown in figure 3(d).

Zoom In Zoom Out Reset image size

We now intend to observe the variation of the MC with respect to the interactions. It is evident from figure 4(a) that the MC predicted from SMF drastically varies from MCS. Utilizing SMF, one can deduce that the MC rises as the strength of interactions is increased which is physically impossible. Clearly, for an interactive system, predictions by SMF largely varies from MCS which can be attributed to the existence of correlation in the system [18] [see figure 4(b)]. Consequently, we need to consider a modified version of SMF that incorporates some of the correlations.

Zoom In Zoom Out Reset image size

To overcome the incapability of SMF for not handling interactions in closed lattice, we employ cluster mean field theory (CMF) which considers some correlations between nearest neighbors. Specifically, we use two-site CMF [2] to factorize the probability of four-site cluster involved in equation (2) to product of probability of two-site cluster as follows:

$$\begin{equation}P\left({\tau }_{i-1},{\tau }_{i},{\tau }_{i+1},{\tau }_{i+2}\right)=P\left({\tau }_{i-1}\vert {\tau }_{i}\right)P\left({\tau }_{i}\vert {\tau }_{i+1}\right)P\left({\tau }_{i+1}\vert {\tau }_{i+2}\right),\end{equation} \tag{ 11 }$$

where

$$\begin{equation}P\left({\tau }_{i}\vert {\tau }_{i+1}\right)=\frac{P\left({\tau }_{i},{\tau }_{i+1}\right)}{P\left({\tau }_{i+1}\right)}.\end{equation} \tag{ 12 }$$

A probability of two-cluster with both empty sites is labeled as P(0, 0) and with two occupied sites is denoted by P(1, 1); whereas half-occupied clusters are labeled as P(1, 0) and P(0, 1). Depending upon the occupancy of sites, each two-cluster can be found in one of the four possible states. Particle–hole symmetry assures that both half-filled two-clusters are equiprobable i.e., P(0, 1) = P(1, 0). Under CMF approximation given by equation (11), the current-density relation is obtained as follows:

$$\begin{equation}J=\frac{\lambda \left(x\right)P\left(1,0\right)}{P\left(1\right)P\left(0\right)}\left[P\left(1,0\right)\left(qP\left(1,0\right)+P\left(0,0\right)+P\left(1,1\right)\right)+rP\left(1,1\right)P\left(0,0\right)\right].\end{equation} \tag{ 13 }$$

By Kolmogorov conditions, we obtain

$$\begin{equation}P\left(1,0\right)+P\left(1,1\right)=\rho ,\end{equation} \tag{ 14 }$$

$$\begin{equation}P\left(1,0\right)+P\left(0,0\right)=1-\rho .\end{equation} \tag{ 15 }$$

Further, from steady state master equation for P(1,0) we procure the following

$$\begin{equation}rP\left(1,1,0,0\right)-qP\left(0,1,0,1\right)=0.\end{equation} \tag{ 16 }$$

Using equations (14)–(16), the probabilities of the two-clusters are expressed as

$$\begin{equation}P\left(1,1\right)=\begin{cases}\rho +\frac{r-\sqrt{r\left(r+4\rho \left(q-r\right)\left(1-\rho \right)\right)}}{2\left(q-r\right)},\quad q,r\ne 1\quad \hfill \\ {\rho }^{2},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad q=r=1\quad \hfill \end{cases}\end{equation} \tag{ 17 }$$

$$\begin{equation}P\left(1,0\right)=\begin{cases}\frac{-r+\sqrt{r\left(r+4\rho \left(q-r\right)\left(1-\rho \right)\right)}}{2\left(q-r\right)},\quad q,r\ne 1\quad \hfill \\ \rho \left(1-\rho \right),\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \enspace q=r=1\quad \hfill \end{cases}\end{equation} \tag{ 18 }$$

$$\begin{equation}P\left(0,0\right)=\begin{cases}1-\rho +\frac{r-\sqrt{r\left(r+4\rho \left(q-r\right)\left(1-\rho \right)\right)}}{2\left(q-r\right)},\quad q,r\ne 1\quad \hfill \\ {\left(1-\rho \right)}^{2},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \enspace q=r=1.\quad \hfill \end{cases}\end{equation} \tag{ 19 }$$

On solving equations (13), (17)–(19), we obtain the following expression for steady state current

$$\begin{align}\hfill J& =\lambda \left(x\right)\left(\frac{\sqrt{{r}^{2}+4r\rho \left(q-r\right)\left(1-\rho \right)}-r}{4{\left(q-r\right)}^{3}\rho \left(\rho -1\right)}\right)\left(4r\rho \left(q-1\right)\left(q-r\right)\left(\rho -1\right)\right.\hfill \\ \hfill & \quad \left.+\left(\sqrt{{r}^{2}+4r\rho \left(q-r\right)\left(1-\rho \right)}-r\right)\left(2rq-q-r\right)\right),\hfill \end{align} \tag{ 20 }$$

which can be replaced in equation (5). Thus, in the steady state, the density profile can be obtained by

$$\begin{align}\hfill {\rho }_{{\pm}}\left(x\right)& =\frac{1}{2}\left[1{\pm}\left(1+\frac{2}{4{r}^{3}\lambda {\left(x\right)}^{2}{\left(q-1\right)}^{2}}\left({J}^{2}{\left(r-q\right)}^{3}+2J\lambda \left(x\right)r\left(q-r\right)\right.\right.\right.\hfill \\ \hfill & \quad \left.\left.{\times}\left(q-3qr+2r\right)-\enspace \lambda {\left(x\right)}^{2}{r}^{2}\left(r-q+2qr\right)+\left(J{\left(r-q\right)}^{2}-\lambda \left(x\right)r\right.\right.\right.\hfill \\ \hfill & \quad \left.{\left.\left.{\times}\left(q+r-2qr\right)\left({J}^{2}{\left(r-q\right)}^{2}-\enspace Jr\lambda {\left(x\right)}^{2}{\left(2q-6r+8qr+r\lambda \left(x\right)\right)}^{\frac{1}{2}}\right)\right)\right)}^{\frac{1}{2}}\right].\hfill \end{align} \tag{ 21 }$$

for q, r ≠ 1. For q = r = 1, we obtain

$$\begin{equation}{\rho }_{{\pm}}\left(x\right)=\frac{1}{2}\left[1{\pm}\sqrt{1-\frac{4J}{\lambda \left(x\right)}}\enspace \right].\end{equation} \tag{ 22 }$$

Here ρ₋(x) ⩽ 0.5 and ρ₊(x) ⩾ 0.5 for all x. The above expressions reduces to the results in reference [25] for λ(x) = 1. We further explore the current for the extreme cases. For E → ∞, P(1, 1) → ρ yields J → 0 which is expected owing to the fact that the particles tend to form clusters in presence of large attractive energy which effectively blocks their movement. In the limit E → −∞, P(1, 1) → 0 which further leads to $J\to \lambda \frac{\rho \left(1-2\rho \right)}{1-\rho }$ (see appendix A). When λ(x) = 1, the expression matches with that of non interacting dimers [28].

Defining X = ρ(1 − ρ) and substituting in equation (20) gives

$$\begin{align}\hfill J& =\lambda \left(x\right)\left(\frac{r-\sqrt{{r}^{2}+4r\left(q-r\right)X}}{4{\left(q-r\right)}^{3}X}\right)\left(4r\left(1-q\right)\left(q-r\right)X\right.\hfill \\ \hfill & \quad \left.+\enspace \left(\sqrt{{r}^{2}+4r\left(q-r\right)X}-r\right)\left(2rq-q-r\right)\right).\hfill \end{align} \tag{ 23 }$$

In order to obtain extrema of current-density relation, we determine ρ such that

$$\begin{equation}\frac{\mathrm{d}J}{\mathrm{d}\rho }=\frac{\mathrm{d}J}{\mathrm{d}X}\left(1-2\rho \right)=0.\end{equation} \tag{ 24 }$$

Clearly, ρ = 0.5 is a critical point for all values of E. Other critical points exist if

$$\begin{equation}{\left(r+\sqrt{{r}^{2}+4r\left(q-r\right)X}\right)}^{2}=\frac{2r\left(q+r-2qr\right)}{1-q}.\end{equation} \tag{ 25 }$$

Substituting X = 0.25 corresponding to ρ = 0.5 in equation (25) yields

$$\begin{equation}{e}^{\theta E}=\frac{1-{e}^{\frac{E}{2}}}{3+{e}^{\frac{E}{2}}}.\end{equation} \tag{ 26 }$$

The above equation can be solved to obtain the critical interaction energy E_c(θ) for 0 < θ ⩽ 1. Note that E_c(θ) remains negative for all 0 < θ ⩽ 1 which means that critical interaction energy is always repulsive in nature. For E < E_c(θ), equation (25) has two real roots. Therefore equation (24) has only one critical point when E ⩾ E_c(θ), and three critical points when E < E_c(θ). For θ = 0.5, equation (26) gives ${E}_{\mathrm{c}}=2\enspace \mathrm{ln}\left(\sqrt{5}-2\right)\approx -2.885{k}_{\mathrm{B}}T$ which coincides with value obtained in [25,29]. Hereafter, we set θ = 0.5 since it splits the interaction energy symmetrically. Depending upon the number of extreme points of equation (20), we further categorize our study into two cases, namely E ⩾ E_c and E < E_c.

2.1.  E ⩾ E_c

For interaction energy greater than or equal to E_c, 0.5 is the only extreme point of current-density relation given by equation (20). Moreover, J attains maximum value at ρ = 0.5 which is given by

$$\begin{equation}{J}_{\mathrm{max}}=\begin{cases}\frac{\lambda \left({x}_{0}\right)\left(\sqrt{qr}-r\right)\left(qr\left(q+r-2\right)+\sqrt{qr}\left(q+r-2qr\right)\right)}{{\left(q-r\right)}^{3}},\quad q,r\ne 1\quad \hfill \\ \frac{\lambda \left({x}_{0}\right)}{4},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad q=r=1,\quad \hfill \end{cases}\end{equation} \tag{ 27 }$$

where λ(x₀) = λ_min. Figure 4(a) shows that the current which is obtained from above expression agrees well with MCS results and overcomes the drawback of SMF approach. Additionally, the density profiles computed using CMF and MCS are in well agreement. (see figure 3).

We now investigate the effect of interactions on the phase diagram in (n − λ_min) plane. It is observed that phase diagram obtained for the proposed model is qualitatively equivalent to that of a non-interacting system [17] with three distinct phases: LD, HD and MC (see figure 5). Also, the density profiles in MC phase exhibits a jump discontinuity and hereafter we refer to it as shock phase (denoted by S_MC). In order to obtain the phase diagram theoretically, we utilized PNC to obtain the phase boundaries. For a fixed λ_min, the boundary between phases LD and S_MC is given by

$$\begin{equation}{n}_{0}={\int }_{0}^{1}{\rho }_{-}\left(x,{J}_{\mathrm{max}}\right)\mathrm{d}x.\end{equation} \tag{ 28 }$$

Similar arguments hold for the phase boundary between S_MC and HD.

Zoom In Zoom Out Reset image size

To comprehend the effect of varying the interaction energy, we first examine the behavior of MC. As E increases, J_max increases until it attains its highest value at a critical energy (${E}_{{\mathrm{c}}_{1}}$) followed by a continuous decrease (see inset, figure 6). For λ(x) = (x − 0.5)² + 0.5, ${E}_{{\mathrm{c}}_{1}}\approx -1.36$ which is theoretically computed from equation (27) and agrees well with that obtained by MCS. Thus, as we increase E, the shock phase shrinks until J_max reaches its maxima, thereafter it starts expanding (see figure 6). This shrinkage and expansion indicates the existence of an interaction energy that has exactly same phase boundaries as that of a non-interacting system. Using equation (28), this interaction energy is obtained to be −2.45 for λ(x) = (x − 0.5)² + 0.25 [see figure 7(a)] and the current in system with E = −2.45 is higher than its non-interactive counterpart. Furthermore, the phase diagram in (n − λ_min) plane for interaction energy E = −2.45 matches well with E = 0 when λ(x) = (x − 0.5)² + λ_min [see figure 7(b)].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2.  E < E_c

When interaction energy is lower than E_c, the current-density relation prescribed by equation (20) admits three distinct extreme points: 0.5, ${\rho }_{{\mathrm{c}}_{1}}$ (<0.5) and ${\rho }_{{\mathrm{c}}_{2}}$ (>0.5), where ${\rho }_{{\mathrm{c}}_{1}}$ and ${\rho }_{{\mathrm{c}}_{2}}$ are the roots of equation (25). Furthermore, J achieves local maximum at ${\rho }_{{\mathrm{c}}_{1}}$ and ${\rho }_{{\mathrm{c}}_{2}}$, and local minimum at 0.5.

With an intend to investigate the effect of interactions, we inspect the phase diagram in n − λ_min plane. Utilizing CMF approach, similar to the case of E > E_c, we obtain three different phases, namely, LD, HD and S_MC. Densities in these phases, obtained by using equation (21), are as follows.

LD phase: the density is smooth throughout the system and is given by ρ(x) = ρ₋(x) where $\rho \left(x\right){\leqslant}{\rho }_{{\mathrm{c}}_{1}}\enspace \forall \enspace x$ with maxima at x₀. At the boundary of LD and S_MC, $\rho \left({x}_{0}\right)={\rho }_{{\mathrm{c}}_{1}}$.

S_MC phase: the density exhibits two shocks, located at x₀ and x_s (>x₀) for which the profile is given by

$$\begin{equation*}\rho \left(x\right)=\begin{cases}_{-}\left(x\right),\quad x{\leqslant}{x}_{0}\quad \hfill \\ {\rho }_{+}\left(x\right),\quad {x}_{0}{< }x{\leqslant}{x}_{s}\quad \hfill \\ {\rho }_{-}\left(x\right),\quad x{ >}{x}_{s}.\quad \hfill \end{cases}\end{equation*}$$

HD phase: the density is smooth throughout the system and is given by ρ(x) = ρ₊(x) where $\rho \left(x\right){\geqslant}{\rho }_{{\mathrm{c}}_{2}}\quad \forall \enspace x$ with minima at x₀. At the boundary of HD and S_MC, $\rho \left({x}_{0}\right)={\rho }_{{\mathrm{c}}_{2}}$.

The theoretically obtained density profiles are validated with MCS for specific set of parameters as shown in figure 8(a). It is evident that in HD and LD phases, the density profiles obtained by CMF and MCS are in good agreement. However, in S_MC phase, MCS reveals one shock which is in contrast to the theoretical findings which predicted the presence of two shocks [see figure 8(b)]. Furthermore, MCS predicts that shock region can be divided into two distinct phases depending upon the position of critical point (ρ = 0.5). In one phase, position of critical point is not fixed, whereas in the other phase, the critical point is fixed at x₀ i.e. ρ(x₀) = 0.5 which has characteristics similar to the S_MC phase obtained when E > E_c. This feature has not been captured by the theoretical finding. We further analyze the shock phase in detail.

Zoom In Zoom Out Reset image size

Breakdown of CMF approach and analysis of shock phase: as discussed above, CMF and MCS results do not agree largely in shock phase. The discrepancy is probably due to the non-homogeneity in density profile that might have increased the correlations which are not captured by proposed CMF theory. This sort of mismatch is also reported in reference [25]. Henceforth we use MCS for analyzing the shock phase.

Now, we inspect the phase transitions for a fixed λ_min by varying n. To analyze the transition from LD to shock phase, we use the following notations for the specific values of n that separate the distinct phases: ${n}_{{b}_{1}}$ and ${n}_{{b}_{2}}$ as shown in figure 9. As n increases from ${n}_{{b}_{1}}$, we observe that the position of ${\rho }_{{\mathrm{c}}_{1}}$ shifts from x₀ to x*(>x₀) and density profile contains only one shock (see figure 10). This is in contrast to findings reported in reference [17] wherein such a shift is not observed. The shifting continues until $n={n}_{{b}_{2}}$ and in that stage, ρ(x₀) attains 0.5 (see figure 11.). It is also noticed that the current in the system also increases with n while ${n}_{{b}_{1}}{< }n{< }{n}_{{b}_{2}}$. This clearly implies that even when the system is in shock phase for ${n}_{{b}_{1}}{< }n{< }{n}_{{b}_{2}}$, it does not attain MC (see figure 12). We denote this phase by S. When n is increased further, the position at which the density profile achieves ${\rho }_{{\mathrm{c}}_{1}}$ and 0.5 become fixed and the current in the system also attains a constant value (see figure 12). This phase is denoted by S_MC. Transition from S_MC to HD through S can be understood in similar lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, depending on current, we observe two types of phases involving a shock:

• (a)  
  When current varies with number of particles while displaying a shock in density profile.
• (b)  
  When current remains constant and density exhibits a shock.

Therefore, in the phase diagram in n − λ_min plane, four distinct phases, namely, LD, HD, S and S_MC are observed. The newly observed phase S has not been reported earlier [17]. Clearly, as opposed to CMF, the MCS current does not attain a minima at n = 0.5 as predicted by CMF (see figure 12).

In previous studies, it has been observed that further decreasing the interaction energy leads to two peaks in current and a local minima at n = 0.5 [25]. On similar lines, we further analyzed the current with respect to n for larger repulsions. Due to computational limitations, we adopted the system length L = 100 and observed that the MCS current attains the double peak with the relative minima at n = 0.5 (see figure 13). Decreasing E further reduces the current, which is expected, owing to fact that as E → −∞, the system converges to that of non-interacting dimers. Formation of double peak can be explained as follows: for dimers, n varies from 0 to 0.5, and the current achieves single peak. By particle–hole symmetry, similar arguments holds for n from 0.5 to 1. Thus, the graph for current versus n has two peaks and is symmetric with respect to n = 0.5.

Zoom In Zoom Out Reset image size

To summarize, we considered a closed TASEP with interacting particles and quenched hopping rates which are further altered depending upon interaction energy E in a thermodynamically consistent manner. We utilized SMF approach which neglects all correlations in the system. It is observed that this theory fails to predict the density profiles and current in the system which is attributed to the correlations that exist in the system. Therefore, to incorporate some correlations, the system is analyzed theoretically using CMF approximation. Specifically, we employed two-site CMF and obtained a critical energy E_c such that the current-density relation has one and three extreme points for E ⩾ E_c and E < E_c, respectively.

It emerges that beyond E_c, results in our proposed model have a behavior qualitatively similar to non-interacting system with three distinct phases: LD, HD and S_MC. Nevertheless, it is interesting to note that as interaction energy increases, S_MC phase reveals a non-monotonic nature i.e, S_MC region decreases followed by an increase after a critical value ${E}_{{\mathrm{c}}_{1}}$. Our observations led us to deduce that there exists an interaction energy for which the boundaries between phases is identical to its non-interactive analogue whereas current is slightly higher in the interactive system.

Below the critical value E_c, CMF approach predicts three phases: LD, HD and S_MC. Density profiles obtained by CMF and MCS agree well in LD and HD phase. However, in shock phase, density profile obtained by CMF predicts presence of two shocks, whereas MCS result exhibits only one shock and divides the shock phase into two disjoint phases: S_MC and S. Thus, MCS reveals four distinct phases: LD, HD, S_MC and S. The striking feature that separates new found S from S_MC phase is that while the current attains its optimal value in S_MC phase, it is not maximal in S. The observed distinguished phase, which appears in our system due to interplay between quenched hopping rates and interactions, has not been reported in earlier relevant studies. For large repulsions, the current exhibits a double peak with respect to n and attains a relative minimum at n = 0.5.

Throughout the study, we have adopted λ(x) = (x − 0.5)² + 0.5 to validate the results with reference [17], however the proposed methodology is generic and can be used even for discontinuous λ(x) (see appendix B). The proposed work is an attempt to understand the NESSs of the closed and quenched driven diffusive system consisting of interacting particles. The two-site CMF approach accurately captures the physical features of the system with interaction energy higher than critical energy E_c. Below E_c, although two-site CMF is in good agreement with MCS in LD and HD phases, it does not seem to work well in the shock region when compared with MCS results. Subsequently, the results for the interaction energy lesser than E_c has been argued with MCS.

We obtain here the expressions for density profiles for a non-interactive system comprising of l-mers i.e., particles of size l. By setting each hole and each l-mer on L sites, as a hole and a monomer, respectively, on another closed lattice with L' sites, we get a mapping L'/L = 1 − (l − 1)ρ_l , where L' = L − (l − 1)N and ρ_l denotes the density of the system with N l-mers and L sites. This leads to the expression of ρ_l as ρ_l = (L'/L)ρ₁, where ρ₁ denotes the density of the system with N monomers on L' sites. Utilizing this mapping gives J_l = (L'/L)J₁. The expressions for steady state current for monomers from [17] leads us to the following expression

$$\begin{equation*}{J}_{l}=\lambda \left(x\right)\frac{{\rho }_{l}\left(x\right)\left(1-l{\rho }_{l}\left(x\right)\right)}{1-\left(l-1\right){\rho }_{l}\left(x\right)},\end{equation*}$$

for steady state current which is similar to that obtained in reference [34].

In our study, the hopping rate function taken into consideration is λ(x) = (x − 0.5)² + λ_min which is a smooth function. It is interesting to note that our analysis is applicable to the system whose hopping rate function might not be differentiable. In fact, our theory also works for a function with discontinuity. However, a shock may appear in the LD and HD phases due to the hopping rate function being discontinuous. We illustrate this in the following subsections.

B.1. Finitely many discontinuities

Figure B1 shows the density profile of a function with two discontinuities. It is evident that system in figure B1(a) has HD phase whereas the system in figure B1(b) has LD phase. Two shocks are present in the density profile. However, these shocks do not imply shock phase. The density profiles show very good agreement of CMF and MCS.

Zoom In Zoom Out Reset image size

B.2. Infinitely many discontinuities

The question that immediately comes to one's mind is: at most how many discontinuities can a function have so that our analysis predict accurate results? As seen in figure B1, it can be safely said that our system works for finite number of discontinuities. But, can we say something about the behavior if the function has infinitely many discontinuities? To understand such a case, consider a hopping rate function λ(x) defined as

$$\begin{equation*}\lambda \left(x\right){:=}\begin{cases}x+0.5,\quad \quad \quad \quad \enspace x=\frac{1}{n},n\in \mathbb{N},\quad \hfill \\ {\left(x-0.5\right)}^{2}+0.5,\quad x\ne \frac{1}{n},n\in \mathbb{N}.\quad \hfill \end{cases}\end{equation*}$$

This function has discontinuities at infinitely many points. We consider a lattice of finite size L(≫1). Then, the hopping rates for each site i is calculated as ${\lambda }_{i}=\lambda \left(\frac{i}{L}\right)$. However, it is worth noting that these λ_i can also be represented as

$$\begin{equation*}\lambda \left(x\right){:=}\begin{cases}x+0.5,\quad \quad \quad \quad \enspace x=\frac{1}{n},\enspace n\in \left\{1,2,\dots ,L-1\right\},\quad \hfill \\ {\left(x-0.5\right)}^{2}+0.5,\quad x\ne \frac{1}{n},n\in \left\{1,2,\dots ,L-1\right\}.\quad \hfill \end{cases}\end{equation*}$$

which again has a finite number of discontinuities and hence, CMF and MCS density profiles match (see figure B2). This is due to the fact that our function with infinite discontinuities reduces to its analogue having finite discontinuities and for such a function, our theory works. In applications, we have finite number of sites which naturally imply that we shall never come across such a λ(x).

Zoom In Zoom Out Reset image size